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Abstract. The BOUT++ code is used to simulate ELMs in a shifted circle 
equilibrium. Reduced ideal MHD simulations are first benchmarked against the 
linear ideal MHD code ELITE, showing good agreement. Diamagnetic drift 
effects arc included finding the expected suppression of high toroidal mode number 
modes. Nonlinear simulations are performed, making the assumption that the 
anomalous kinematic electron viscosity is comparable to the anomalous electron 
thermal diffusivity. This allows simulations with realistically high Lundquist numbers 
(S = 10 8 ), finding ELM sizes of 5-10% of the pedestal stored thermal energy. Scans 
show a strong dependence of the ELM size resistivity at low Lundquist numbers, 
with higher resistivity leading to more violent eruptions. At high Lundquist numbers 
relevant to high-performance discharges, ELM size is independent of resistivity as 
hyper-resistivity becomes the dominant dissipative effect. 



PACS numbers: 52.25.Xz, 52.65.Kj, 52.55.Fa 
1. Introduction 

Future tokamak devices such as ITER and DEMO require high performance plasma 
operation in order to demonstrate economical fusion performance. For this reason, 
discharges with a transport barrier close to the plasma edge (H-mode pQ) are the current 
baseline operating scenario for ITER [2]. Whilst transport barriers at the plasma edge 
result in improved performance, the steep pressure gradients and associated bootstrap 
current can destabilise peeling-ballooning modes [3J. These ideal MHD instabilities are 
thought to be responsible for triggering the observed eruptions of filamentary structures 
from the plasma edge known as Edge Localised Modes (ELMs). The particles and energy 
released during ELMs are deposited on material surfaces and are potentially damaging 
in future devices. There is therefore interest in understanding and controlling these 
events. 

Understanding of the linear onset and structure of peeling-ballooning modes is now 
quite well developed, with codes such as ELITE $ H], GATO [5j and MISHKA [7] 
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able to predict experimental operating limits[3l El E]. Study of the nonlinear evolution 
of these instabilities is more challenging, and although there are analytic theories [10] 
and semi-analytic models [HJ [12] for the early non-linear evolution, it is not yet fully 
understood how this develops into the filamentary structures observed, and ultimately 
how particles and energy are lost. 

Several 3D non-linear codes have been used to simulate ELMs, including NIMROD 
[131 [HI HS1 US] , BOUT [HI CEE], JOREK OEniEI], GEM [221123], M3D [2H23I2S] and 
M3D-C1 [27]. These codes incorporate a wide range of physics including (in the case 
of BOUT and some NIMROD simulations) 2-fluid effects. In this paper, simulations of 
ELMs using the the BOUT++ code [2H] based on modified reduced MHD in a shifted 
circle equilibrium are presented, expanding on and extending results presented in [29J. 
By introducing a hyper-resistive term into Ohm's law, realistic resistivities have been 
simulated, and the effect of varying resistivity on ELM sizes has been investigated. 

Section [2] describes the starting set of equations used for comparison with linear 
ideal MHD, and the equilibrium simulated is described in section [3j A novel means to 
handle the vacuum region which has been found to work well in ideal or resistive MHD 



simulations is described in section 2.1 Linear simulations are benchmarked against 
linear ideal MHD codes in section [4], and nonlinear simulations of ideal reduced MHD 
and the issues encountered are discussed in section [5j By incorporating diamagnetic 
drift and either high resistivity or a hyper-resistivity, simulations of nonlinear eruptions 
are performed and discussed in sections [6] and [JJ respectively. The effect of varying 
resistivity on ELM crash size is then discussed in more detail in section [8j Details 
of the magnetic field structure are presented in section [9j indicating a reconnection of 
flux-surfaces at the leading edges of the erupting filaments. 



2. Model equations 

In this paper, the starting equations are those of high-/? reduced MHD [30], evolving 
pressure P, vorticity U and magnetic potential if) = A\\/Bq\ 
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where V\\F = Bd\\ (F/B) for any F, d\\ = (b + bi) • V, b = B/B. The plasma density 
is constant in space and time, fixed at 10 19 m -3 . The perturbed magnetic field is given 
by B x = Vip x B where B is the unperturbed field. k = b • Vb is the field curvature 



Simulation of ELMs using BOUT++ 



3 



and all quantities are in SI units. The vorticity equation [T] includes the kink/peeling 
term through the perturbed magnetic field: 



This minimal set of equations contains the basic physics needed to describe peeling- 
ballooning modes, including pressure/curvature and parallel current (kink/peeling) 
instability drives, and field-line bending stabilisation. The model is a useful starting 
point because it allows benchmarking against linear ideal MHD codes such as ELITE, 
providing a valuable means of checking the results of these simulations. 

In the pedestal with steep pressure gradients diamagnetic effects are expected to 
be significant, and damp high toroidal mode-number, n, instabilities which can lead to 
problems in ideal MHD simulations (see section [5]). Therefore, the set of equations (JTjjs]) 
is modified in section [6] to include the diamagnetic drift. 

Dirichlet (zero value) boundary conditions are used for the perturbed pressure, 
vorticity and parallel current. To be consistent, the boundary conditions for ip and <p 
must satisfy V^V = an d = 0. This is done by solving for each toroidal mode- 
number which have solutions n ~ C+exp C^fg^Jg^k^x) + C_ exp (—yg^Jg^k^ 
where ( is the toroidal angle, = n/R and x is the radial (ip) coordinate. Only 
solutions which decay going out of the domain are allowed, so C + = on the outer 
boundary, and C_ = on the inner boundary. 

2.1. Vacuum region 

A complication of simulating the plasma edge region is the handling of the vacuum 
region, and the plasma-vacuum interface. In linear codes such as ELITE an analytical 
approach can be used; non-linear codes must also handle motion of the vacuum interface. 
The usual approach has been to treat the vacuum as a resistive plasma, with a jump in 
resistivity between plasma and vacuum [16]. Here a different approach is used: in the 
vacuum we evolve ip towards a self-consistent solution determined by currents in the 
core only. The procedure is as follows. 

First define a smoothed step function 0, which switches between in the core to 1 
in the vacuum: 



where P vac is the pressure at the plasma- vacuum interface, and AP vac < P vac is the 



which may or may not include currents in the vacuum region. From this, a "target" 
current j^ ar ' 9et [ s calculated by setting all vacuum currents to zero 



b- V 
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This is the closest physically acceptable solution to the current result. From this, a 
target if) can be calculated 

garget = y-2 fajtW* / Bq ) 

In the vacuum region, this is the solution which would give zero current, but simply 
setting if) to this value in the vacuum would result in an inconsistency because if> in the 
core is affected by currents in the vacuum. Instead, equation [2] is modified to 

^ = - (1 " ©) ^V|,0 + e - ip) /r jvac 

and so if) converges on the target value with a small time-constant Tj vac . if) therefore 
evolves to a self-consistent state with zero current in the vacuum region. 

This method has been found to work well for ideal and resistive simulations, even 
in the nonlinear regime. It is useful because this method uses the same code already 
used for solving the reduced MHD model with minimal modifications. At this point this 
method does not work once diamagnetic effects are included and so is not used in the 
nonlinear results presented in sections [6] and [7j 



3. Equilibrium 

Simulation of ELMs in full x-point geometry is necessary in order to predict the 
behaviour of ELMs in future devices. As a starting point however, it is useful to use 
equilibria which remove many of the complications of real magnetic geometry: This 
simplifies study of the basic physics of ELMs, but more importantly enables accurate 
benchmarking of the results against linear theory and other codes. 

A test case which is becoming standard in this field is the cbml8 series, created 
by P.Snyder using the TOQ equilibrium code. The case used here cbml8_dens8 has 
now been used by NIMROD |i6j and M3D-cl [27] to benchmark linear growth-rates 
against the ELITE [3l H] and GATO [5] linear MHD stability codes. The pressure and 



current profiles are shown in figure 1 (a) over the range of normalised if) simulated using 
BOUT++ (0.4 < ip n < 1.2). The definition of the plasma edge is somewhat arbitrary 
since here there are no separatrices, and so ip n — 1 is defined as the point where the 
equilibrium plasma pressure gradient falls to zero. For this equilibrium, the normalised 
beta is (3 N = 1.51, edge q a = 3.03, pedestal toroidal pressure f3 t0 = 1.9 x 10 -2 and 
normalised pedestal width is L pe d/a = 0.049. 

The mesh used for the simulations in this paper is shown in figure l(b)| This has 



a resolution of 512 points in radial coordinate if> and 64 in poloidal coordinate 9. In the 
toroidal direction 64 points are usually used, with the results checked by doubling this 
to 128. Note that despite this relatively low resolution in 9, high m numbers can be 
simulated. This is because the field-aligned coordinate system 

,0 

x = ip — if>o y = 9 z = cj) — I v (if), 9) d9 

Je 

(where v{if),9) is the field-line pitch) means that the 9 direction is aligned with the 
magnetic field and so hq = 64 is therefore the number of points along a field-line 
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Figure 1. Equilibrium (cbml8_dens8) profiles and mesh 



per poloidal turn. The effective poloidal resolution depends on the toroidal resolution 
(number of field-lines simulated x number of repetitions in a full torus) and the local 
safety factor q ~ q a — 3.03, and so is ~ 3 x 64 x 5 = 960 for the simulations shown here. 

In order to avoid problems associated with magnetic shear and deformation of 
grid cells associated with field-aligned coordinate systems, shifted local coordinate 
systems [211 122] are used so that the radial derivative is always taken in ip. Further 
details of this coordinate system and its implementation in BOUT++ can be found in 



4. Linear benchmarking 

In order to benchmark BOUT++, linear simulations have been performed and 
comparison made to the ELITE [21 E] and GATO [5j linear ideal MHD codes. Figure [2] 
shows the linear growth-rate from BOUT++ for this equilibrium (star symbols), along 
with the results from ELITE (open circles) and GATO (filled squares). All growth-rates 
are normalised to an Alfven frequency ua = Va/R , with Va the Alfven velocity, and R 
the major radius. For all the simulations presented in this paper, a reference density of 
10 19 m~ 3 was used, so that time is normalised to ta = 0.37/is. 

At high-n, the BOUT++ result begins to deviate from the ELITE result, whereas 
it might be expected that finite-n corrections would lead to the greatest differences 
being at low n. The reason for this is that as the mode number n is increased the 
radial width of the mode becomes narrower: Since these simulations all used the same 
radial resolution, at high n small-scale structures can no longer be resolved properly. 
The mode structure of the radial displacement £^ from ELITE and BOUT++ have 
also been compared [28], finding good agreement. Individual (m, n) modes are found to 
peak close to their resonant magnetic surfaces, as is expected from analytic theory. This 
test is a proof-of-principle which demonstrates that BOUT+-I- is capable of simulating 
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Figure 2. Linear growth rate of a strongly ballooning test calculated by 

BOUT++ and the linear ideal MHD codes ELITE and GATO 



the ideal ballooning mode correctly using ideal reduced MHD. Linear benchmarks are a 
necessary but not sufficient condition for non-linear results: comparisons against other 
nonlinear codes and ultimately experimental results are needed to verify the fidelity of 
nonlinear results. 

5. Non-linear ideal MHD 

In ideal MHD linear simulations, peaks in the parallel current form at rational q surfaces. 
The width of these current sheets is of the order of the distance between rational surfaces 
A q ~ r/nqs ~ 6mm, which can be resolved in the linear regime since when using 512 
radial points across the pedestal the grid spacing is ~ 1.1mm. Linear codes such as 
ELITE and MARS [33] typically pack mesh-points around rational surfaces to handle 
these features, but this becomes much more challenging in the nonlinear regime: flux- 
surfaces move over time, current sheets can be compressed so reducing scale lengths, 
and any adaptive scheme will struggle to handle highly nonlinear regimes where flux 
surfaces may be destroyed. 

Nonlinear simulations of ideal MHD have been attempted using BOUT++, but as 
yet without success: As flux-surfaces are distorted, small-scale structures are generated 
in the binormal direction which corresponds to the generation of high toroidal mode- 
numbers. The size of the structures formed in the toroidal direction is approximately 
given by rotating the radial scales in the perpendicular plane, and projecting in the 
toroidal direction which gives ~ AB^/Bg ~ 5cm and so a maximum toroidal mode- 
number of n m ax ~ 600. 

A second source of high- A; structures and hence numerical problems is that the 
modes with the highest growth-rate are those with the highest n (i.e. those at the grid 
scale). Any perturbation at high n will therefore rapidly grow and eventually dominate 
the solution. 
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Handling small-scale structures is a problem in many nonlinear simulations, and 
if possible then these should be resolved. In practice however, it is often necessary to 
add some sort of numerical dissipation to remove structures at the grid scales. Several 
approaches are commonly used to achieve this. For example the solution can be filtered 
in Fourier space (so only the low-n modes are simulated) or a hyper-diffusion like term 
can be added. Whichever method is used, it should be chosen to minimise the impact 
on the large-scale structures of interest. 

Tokamak simulations have the advantage that Larmor radius effects should 
naturally damp high-fc structures. Adding effects such as diamagnetic drifts both 
make the simulation more realistic and also aid numerical stability by damping high- 
k structures. In this paper, nonlinear simulations are performed by first adding 
diamagnetic drift and following the common practice of incorporating a high resistivity 
in section [6j before introducing a hyper-resistive term in section [7} 



6. Diamagnetic and resistive simulations 

For these initial nonlinear simulations, the basic effects of diamagnetic stabilisation are 
incorporated into the model by modifying the vorticity Q so that the total plasma 
velocity is now given by the sum of E x B and diamagnetic drifts. 

V ~ V E + V D = i-B x V0 + B x VP 

B z enB z 

Hence assuming constant mass density n, this gives 

1^2 f. ■ P\ 



U=oV3U+- (7) 
B \ en) 

In this simplified model, both equilibrium and turbulent zonal (m = n = 0) flows 
have been set to zero: Veo + V^o = and toroidally-averaged (5ve + Svd) = 0. 
The equilibrium radial electric field is therefore given by E r0 = (l/n Zje) V r Pjo (with 
Pio = Pq/2) and perturbed radial electric field (E r ) = (l/no^e) V r (Pi)- Recent work 
[20} I2TI [33] indicates that the nonlinear formation of zonal flows results in multiple 
filaments breaking off from the plasma, but investigation of this effect is beyond the 
scope of this paper and the subject of future work. 

At realistic resistivities of S ~ 10 8 , the scale of the current sheets is limited by 
resistivity to approximately Ar ~ R^J (ujaH) / S ~ 10 — 100/xm, of the same order as 
the electron gyro-radius. The radial grid spacing in these simulations ~ 1mm, and 
so to limit the smallest scales to these sizes a greatly increased resistivity of S = 10 4 is 
used, constant in space. 

Diamagnetic drift and high resistivity produces the linear growth-rates shown in 
figure [3j High resistivity results in an increase in growth-rate above the ideal MHD 
case (open circles), whilst high-n modes are suppressed by diamagnetic drift effects in 
line with theoretical expectations. The peak growth-rate now occurs at n — 15, and so 
this toroidal mode-number is used as a starting point for the nonlinear results shown 
in figure |4j Figure 4(a) shows the plasma displacement, measured by finding where 
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Figure 3. Linear growth rate for simulations with and S = 10 4 
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(c) Plasma energy loss AW pe d/W ped 
Figure 4. Nonlinear simulation results with and S = 10 4 



the pressure crosses a threshold. Outwards displacement is shown as a solid line, and 
inwards displacement as a dashed line: These are approximately equal in the linear 
regime (as expected since mode goes like exp (in() ), but begin to diverge once the 
displacement exceeds ~ 1 — 2cm with the outgoing filament growth-rate reducing whilst 
the ingoing hole continues to grow at the linear rate. 
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Figure 4(b) shows the pressure profiles averaged over (equilibrium) flux surfaces 



at the times indicated by vertical lines in figure 4(a) . This shows that this eruption 



propagates far into the plasma beyond the steep pressure-gradient region (shown in the 



t = equilibrium profile in figure 4(b) , and in figure 1(a) ). 

To quantify the size of the ELM eruption, and provide a measure which can be 
compared against typical experimentally measured ELM sizes, the thermal energy from 
the inner edge of the domain (R = 4.2m at the outboard mid-plane) to R = 4.6m 
(vertical line in figure [4(b) ) W pe a = J 3P/2dV is used. The fraction of this energy lost 



during the ELM crash AW pe d/W pe d is shown in figure 4(c) For this highly resistive 
case, approximately 30% of the thermal energy in the pedestal is lost, limited by the 
radial size of the computational domain. 

7. Diamagnetic and hyper-resistive simulations 

Using high resistivities to damp small-scale structures and maintain numerical stability 
in the nonlinear regime results in greatly increased growth-rates, modified stability 
thresholds, and large eruptions which propagate far into the plasma. 

An alternative effect which damps high-fc structures whilst minimising the effect on 
large-scale structures is hyper-resistivity [291 H23 [36] : The parallel electric field equation 
for if) = An /Bo is modified to 

at .dq /i /io 

where rjH is a hyper-resistivity. The hyper-Lundquist number is then defined analogously 
to the Lundquist number as the dimensionless ratio of an Alfven wave crossing time 
to a hyper-diffusion timescale. Sh = HoRoVa/vh = S/och- The hyper-Lundquist 
parameter an = f]H/R 2 ^ei can be estimated from electron collisional viscosity [2H] as 
an — ^e/R 2 v e i if it is assumed that the electron viscosity /i e is comparable to the 
anomalous electron thermal diffusivity fi e ~ \ e ~ lm 2 /s. Taking a value of v ei ~ 10 5 s _1 
gives an — 10~ 4 — 1CT 6 for R ~ 3m. 

Taking S = 10 8 and Sh = 10 12 , the perpendicular scale-length of the parallel 
current can be estimated this time by balancing dip/dt against the hyper- resistive term. 
This gives Ah — R (oja/'j/Sh) 1 ^ — 2mm, which is of the order of the radial grid spacing 

~ 1.1mm, and an order of magnitude larger than was the case with only resistivity 
included. By including this hyper-resistive effect, realistic resistivities of S ~ 10 8 can 
be simulated. 

Linear growth-rates are shown in figure [5j with the ideal MHD and resistive 
(figure [3]) results for comparison. Growth-rates are below the ideal MHD case, stabilised 
at high-n by diamagnetic drifts as expected, and recover the low-n ideal MHD growth- 
rates. This is encouraging and an improvement on the resistive case where enhanced 
growth-rates were found. 

Nonlinear simulation results are shown in figure k3] using S = 10 8 and S = 10 12 , 



with the same plots as in figure m As with the resistive case, figure 6(a) shows that 
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Figure 5. Linear growth rate for simulations with u;* and Sh = 10 
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Figure 6. Nonlinear simulation results with a;*, S = 10 s and Sjj — 10 12 



the growth-rate approximately follows the linear growth-rate until a displacement of 
approximately 1 — 2cm (up to t ~ 75r^) at which point the growth-rate of the plasma 
displacement slows and the profiles relax. 

A striking difference between this low-resistivity case and the resistive case 
presented in section [6] is that the eruption no longer penetrates far into the core. 
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Instead, after a fast initial stage the profiles relax on timescales of 10's of ta as shown 
in figure |6(b) where the surface-averages pressure profiles are plotted for the times 



indicated as vertical lines in figure 6 (a" 



The ELM loss AW pe d/W pe d is plotted as a function of time in figure [6] which shows 
a loss of ~ 6% of the pedestal thermal energy (with a continuing slow relaxation of the 
profiles). Most of this is due to convective losses in the initial fast stage of the eruption 
as filaments are ejected from the plasma. Note that the simplified model used here 
does not include parallel heat conduction, and so only convective losses are captured. 
Convective losses are thought to be important in x-point simulations where open field- 
lines can quickly transport away heat, but will probably not play such a key role in 
these circular "limiter" plasmas without open field-lines. 



8. Effect of resistivity on ELM size 

The size of the ELM eruption varies dramatically between the high and low resistivity 
cases (figure [5] and [6] respectively). The transition between these cases is shown in 
figure [7] which plots the ELM loss fraction AW pe d/W pe d as a function of the Lundquist 
number for a fixed Sh = 10 12 . As the Lundquist number is increased, the loss of thermal 
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Figure 7. ELM size against Lundquist number for fixed Sh = 10 12 

energy during the ELM (open circles in figure [7]) is seen to drop from ~ 30% to ~ 6% 
over a range of S = 10 5 — 10 7 . One explanation for the more violent eruptions at high 
resistivity is the increased linear growth-rate: Also plotted in figure [7] (crosses) is the 
linear growth-rate of the n = 15 mode, which also shows an increase over a similar range 
of Lundquist number. 

Since there is little change in the linear growth-rate between S = 10 6 and 5* = 10 7 , 
whilst ELM size changes by a factor of ~ 2, the process which determines ELM loss may 
be a nonlinear phenomenon related to dissipation at the smallest scales. The point at 
which the resistivity starts to dominate over the hyper-resistivity (A# = Ah) is given 
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by S = ^Sh^a/i ~ 2 x 10 6 . Increasing the Lundquist number above 10 6 therefore 
has little effect on the ELM size AW pe d/W pe d because hyper-resistivity is the dominant 
dissipative effect. This may then imply that in high Lundquist number regimes relevant 
to high-performance discharges, the (convective) ELM loss is determined by the hyper- 
resistive dissipation rather than the resistivity. 



9. Magnetic field structure 



In order to understand the ELM eruption in the low resistivity case (figure [6]), Poincare 
(puncture) plots of the magnetic field structure have been produced and compared 
against the pressure contours. Figure [8] shows the (normalised) pressure 2n p/B 2 at the 
outboard mid-plane as a function of if) and toroidal angle ( for t = 60ta and t = 70ta- 
On top of these are plotted a puncture plot of the magnetic field showing the field 
structure at these two times. Up to t ~ 60 — 65r^, the magnetic field surfaces follow 
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Figure 8. Pressure and magnetic field puncture plot at outboard mid-plane. Minor 
radius increases going up the figures 



pressure isosurfaces (frozen-in condition), and the growth-rate remains approximately 
the same as the linear mode (figure 6(a)). During this period, the leading edge of 
the erupting filaments steepens. This forces flux-surfaces closer together at the tips of 
the filaments, reducing perpendicular scale lengths and increasing the importance of 
(hyper-)resistivity. Eventually, the flux-surfaces break and the magnetic field becomes 
disordered. After this time (t ~ 65ta), the eruption growth-rate falls and the profiles 
begin to relax on a slower timescale. 



10. Conclusions 



Simulations of Edge Localised Modes (ELMs) have been performed using the BOUT++ 
simulation code [28J for a shifted circle equilibrium (cbml8_dens8). The linear structure 
and growth-rates of the unstable ballooning modes have been benchmarked without 
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diamagnetic effects against the linear ideal MHD codes ELITE and GATO showing 
good agreement. By incorporating diamagnetic drift and hyper- resistivity, the scale 
length of J|| structures can be resolved, allowing nonlinear simulations with resistivities 
relevant to high-performance H-mode tokamak edge plasmas. These show reasonable 
ELM losses of ~ 6% of the pedestal thermal energy, with a limited radial extent, in 
contrast to simulations with enhanced resistivities. 

Results at high Lundquist number indicate that eruptions consist of two stages: a 
fast emergence of plasma "fingers" at approximately the linear growth-rate (figure 6(a) ), 
during which time the frozen-in condition is obeyed (figure 8(a) ). This eruption results 
in a steepening of the pressure profile at the tips of the plasma fingers, forcing flux 
surfaces closer together and steepening local pressure gradients. Once the fingers have 
emerged ~ 1 — 2cm, reconnection occurs at these high-gradient regions, the growth-rate 
is reduced and filaments begin to break up as the plasma profiles relax (figure 6(b) ). 

Experimentally, it is found that the ELM size AW pe d/W pe( i decreases as collisionality 
v* is increased j2j |37] . The results presented in figure [7] indicate the opposite trend at 
low Lundquist numbers, with the size of ELM eruptions increasing with collisionality 
given the same starting pressure and J\\ profiles. At high Lundquist numbers (i.e. low 
u*) relevant to H-mode pedestals, ELM size has been found to be quite insensitive to 
Lundquist number, with dissipation dominated by hyper-resistivity. The variation of 
hyper-resistivity with collisionality is however not well known. Assuming a constant 
electron viscosity /i e ~ lm 2 /s gives an oc 1/V, and so Sh ~ const (since S oc l/u). This 
may imply that convective ELM losses are approximately constant with collisionality, 
and that the observed variation of ELM loss with v* is due to another process such as 
increasing parallel conductive losses at low collisionality. 

A important consideration when interpreting the results presented here is that in 
these simulations, variation in H-mode pedestal characteristics or bootstrap current 
with collisionality has not been included, though in practice all these things would vary 
and change the point at which an ELM is triggered. Predictions of ELM size in future 
devices depends on a model of the H-mode pedestal from which to start a nonlinear 
ELM simulation, either coupling kinetic and fluid models such as work coupling XGCO 
to M3D [38J. Other possibilities include using a semi- heuristic model such as EPED1 [39J 
to provide input to a starting equilibrium. 

Further work includes the incorporation of parallel heat conduction effects in order 
to study the conductive losses in conjunction with the convective losses studied here. 
As mentioned in section |6j non-linear generation of poloidal flows have been found to be 
important in breaking off ELM filaments [201 EE] and so incorporation of these effects 
into BOUT++ simulations is a priority. 
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